

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>agpy.collapse_gaussfit &mdash; agpy 0.1.2 documentation</title>
    
    <link rel="stylesheet" href="../../_static/extra.css" type="text/css" />
    <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    '../../',
        VERSION:     '0.1.2',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="../../_static/jquery.js"></script>
    <script type="text/javascript" src="../../_static/underscore.js"></script>
    <script type="text/javascript" src="../../_static/doctools.js"></script>
    <link rel="top" title="agpy 0.1.2 documentation" href="../../index.html" />
    <link rel="up" title="agpy" href="../agpy.html" />
     
    <script type="text/javascript">

      var _gaq = _gaq || [];
      _gaq.push(['_setDomainName', 'pyspeckit.bitbucket.org']);
      _gaq.push(['_setAllowHash', false]);
      _gaq.push(['_trackPageview']);


    </script>
    <link rel="stylesheet" type="text/css" href="../../_static/extra.css" />
  </head>
  <body>
    <div class="header-wrapper">
      <div class="header">
        <h1><a href="../../index.html">agpy 0.1.2 documentation</a></h1>
        <div class="rel">
          <a href="http://agpy.googlecode.com">agpy Home </a>  |
          <a href=../../index.html>Docs Home </a>  |
          <a href="http://code.google.com/p/agpy/w/list">Wiki</a>  |
          <a href=../../search.html>Search </a>
        </div>
       </div>
    </div>

    <div class="content-wrapper">
      <div class="content">
        <div class="document">
            
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <h1>Source code for agpy.collapse_gaussfit</h1><div class="highlight"><pre>
<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">-----------------</span>
<span class="sd">Collapse Gaussfit</span>
<span class="sd">-----------------</span>

<span class="sd">This was an early attempt to automate gaussian fitting over a data cube using</span>
<span class="sd">(multiple) gaussian decomposition for each spectrum.   It&#39;s reasonably</span>
<span class="sd">effective, but the uses are somewhat minimal.  I&#39;ve tried shifting my</span>
<span class="sd">cube-related work to `pyspeckit &lt;pyspeckit.bitbucket.org&gt;`_.</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="k">try</span><span class="p">:</span>
    <span class="kn">import</span> <span class="nn">scipy</span>
    <span class="kn">from</span> <span class="nn">scipy</span> <span class="kn">import</span> <span class="n">optimize</span><span class="p">,</span><span class="n">sqrt</span>
    <span class="kn">from</span> <span class="nn">scipy.optimize</span> <span class="kn">import</span> <span class="n">leastsq</span>
    <span class="c">#from scipy.stats.stats import nanmedian,nanmean,_nanmedian</span>
<span class="k">except</span> <span class="ne">ImportError</span><span class="p">:</span>
    <span class="k">print</span> <span class="s">&quot;Scipy cold not be loaded.  Collapse_gaussfit may fail&quot;</span>
<span class="kn">import</span> <span class="nn">numpy</span>
<span class="kn">from</span> <span class="nn">numpy</span> <span class="kn">import</span> <span class="n">vectorize</span><span class="p">,</span><span class="n">zeros</span><span class="p">,</span><span class="n">exp</span><span class="p">,</span><span class="n">median</span><span class="p">,</span><span class="n">where</span><span class="p">,</span><span class="n">asarray</span><span class="p">,</span><span class="n">array</span><span class="p">,</span><span class="n">nonzero</span><span class="p">,</span><span class="n">ma</span><span class="p">,</span><span class="n">arange</span><span class="p">,</span><span class="n">square</span>
<span class="kn">import</span> <span class="nn">matplotlib</span>
<span class="c">#matplotlib.use(&#39;Agg&#39;)</span>
<span class="kn">from</span> <span class="nn">pylab</span> <span class="kn">import</span> <span class="n">indices</span><span class="p">,</span><span class="n">figure</span><span class="p">,</span><span class="n">clf</span><span class="p">,</span><span class="n">savefig</span><span class="p">,</span><span class="n">plot</span><span class="p">,</span><span class="n">legend</span><span class="p">,</span><span class="n">text</span><span class="p">,</span><span class="n">axes</span><span class="p">,</span><span class="n">title</span>
<span class="kn">import</span> <span class="nn">pickle</span>
<span class="kn">import</span> <span class="nn">pyfits</span>
<span class="kn">import</span> <span class="nn">time</span>
<span class="kn">from</span> <span class="nn">mad</span> <span class="kn">import</span> <span class="n">MAD</span>
<span class="kn">from</span> <span class="nn">ratosexagesimal</span> <span class="kn">import</span> <span class="n">ratos</span><span class="p">,</span><span class="n">dectos</span>

<div class="viewcode-block" id="nanmedian"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.nanmedian">[docs]</a><span class="k">def</span> <span class="nf">nanmedian</span><span class="p">(</span><span class="n">arr</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot; nanmedian - this version is NOT capable of broadcasting (operating along axes) &quot;&quot;&quot;</span>
    <span class="k">return</span> <span class="n">median</span><span class="p">(</span><span class="n">arr</span><span class="p">[</span><span class="n">arr</span><span class="o">==</span><span class="n">arr</span><span class="p">])</span></div>
<div class="viewcode-block" id="nanmean"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.nanmean">[docs]</a><span class="k">def</span> <span class="nf">nanmean</span><span class="p">(</span><span class="n">arr</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot; nanmean - this version is NOT capable of broadcasting (operating along axes) &quot;&quot;&quot;</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">arr</span><span class="p">[</span><span class="n">arr</span><span class="o">==</span><span class="n">arr</span><span class="p">])</span><span class="o">.</span><span class="n">mean</span><span class="p">()</span>

<span class="c"># read in file</span>
<span class="c"># filename = sys.argv[1]</span>
<span class="c"># fitsfile = pyfits.open(filename)</span>
<span class="c"># cube = fitsfile[0].data</span>

<span class="c"># def gaussian(dx,sigma):</span>
<span class="c">#     return lambda x: exp( - (x-dx)**2 / sigma**2 )</span>
<span class="c"># def return_param(xarr,param):</span>
<span class="c">#     errorfunction = lambda p:gaussian(*p)(*indices(xarr.shape))-xarr</span>
<span class="c">#     pars, cov, infodict, errmsg, success = optimize.leastsq(errorfunction, [len(xarr)/2.,1], full_output=1)</span>
<span class="c">#     print errmsg</span>
<span class="c">#     if param == &#39;width&#39;:</span>
<span class="c">#         return pars[1]</span>
<span class="c">#     elif param == &#39;center&#39;:</span>
<span class="c">#         return pars[0]</span>
<span class="c">#     else:</span>
<span class="c">#         return </span></div>
<div class="viewcode-block" id="gaussian"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.gaussian">[docs]</a><span class="k">def</span> <span class="nf">gaussian</span><span class="p">(</span><span class="n">dx</span><span class="p">,</span><span class="n">sigma</span><span class="p">,</span><span class="n">a</span><span class="p">):</span>
    <span class="k">return</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">a</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span> <span class="o">-</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">dx</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">sigma</span><span class="o">**</span><span class="mi">2</span> <span class="p">)</span></div>
<div class="viewcode-block" id="double_gaussian"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.double_gaussian">[docs]</a><span class="k">def</span> <span class="nf">double_gaussian</span><span class="p">(</span><span class="n">dx1</span><span class="p">,</span><span class="n">dx2</span><span class="p">,</span><span class="n">sigma1</span><span class="p">,</span><span class="n">sigma2</span><span class="p">,</span><span class="n">a1</span><span class="p">,</span><span class="n">a2</span><span class="p">):</span>
    <span class="k">return</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">a1</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span> <span class="o">-</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">dx1</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">sigma1</span><span class="o">**</span><span class="mi">2</span> <span class="p">)</span> <span class="o">+</span>  <span class="n">a2</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span> <span class="o">-</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">dx2</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">sigma2</span><span class="o">**</span><span class="mi">2</span> <span class="p">)</span></div>
<div class="viewcode-block" id="triple_gaussian"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.triple_gaussian">[docs]</a><span class="k">def</span> <span class="nf">triple_gaussian</span><span class="p">(</span><span class="n">dx1</span><span class="p">,</span><span class="n">dx2</span><span class="p">,</span><span class="n">dx3</span><span class="p">,</span><span class="n">sigma1</span><span class="p">,</span><span class="n">sigma2</span><span class="p">,</span><span class="n">sigma3</span><span class="p">,</span><span class="n">a1</span><span class="p">,</span><span class="n">a2</span><span class="p">,</span><span class="n">a3</span><span class="p">):</span>
    <span class="k">return</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="nb">abs</span><span class="p">(</span><span class="n">a1</span><span class="p">)</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span> <span class="o">-</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">dx1</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">sigma1</span><span class="o">**</span><span class="mi">2</span> <span class="p">)</span> <span class="o">+</span>  <span class="nb">abs</span><span class="p">(</span><span class="n">a2</span><span class="p">)</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span> <span class="o">-</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">dx2</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">sigma2</span><span class="o">**</span><span class="mi">2</span> <span class="p">)</span> <span class="o">+</span>  <span class="nb">abs</span><span class="p">(</span><span class="n">a3</span><span class="p">)</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span> <span class="o">-</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">dx3</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">sigma3</span><span class="o">**</span><span class="mi">2</span> <span class="p">)</span></div>
<div class="viewcode-block" id="n_gaussian"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.n_gaussian">[docs]</a><span class="k">def</span> <span class="nf">n_gaussian</span><span class="p">(</span><span class="n">dx</span><span class="p">,</span><span class="n">sigma</span><span class="p">,</span><span class="n">a</span><span class="p">):</span>
    <span class="k">def</span> <span class="nf">g</span><span class="p">(</span><span class="n">x</span><span class="p">):</span>
        <span class="n">v</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">dx</span><span class="p">)):</span>
            <span class="n">v</span> <span class="o">+=</span> <span class="n">a</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">*</span> <span class="n">exp</span><span class="p">(</span> <span class="o">-</span> <span class="p">(</span> <span class="n">x</span> <span class="o">-</span> <span class="n">dx</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="n">sigma</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span> <span class="p">)</span>
        <span class="k">return</span> <span class="n">v</span>
    <span class="k">return</span> <span class="n">g</span></div>
<div class="viewcode-block" id="gerr"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.gerr">[docs]</a><span class="k">def</span> <span class="nf">gerr</span><span class="p">(</span><span class="n">xarr</span><span class="p">):</span>
    <span class="k">return</span> <span class="k">lambda</span> <span class="n">p</span><span class="p">:</span><span class="n">xarr</span><span class="o">-</span><span class="n">gaussian</span><span class="p">(</span><span class="o">*</span><span class="n">p</span><span class="p">)(</span><span class="o">*</span><span class="n">indices</span><span class="p">(</span><span class="n">xarr</span><span class="o">.</span><span class="n">shape</span><span class="p">))</span></div>
<div class="viewcode-block" id="double_gerr"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.double_gerr">[docs]</a><span class="k">def</span> <span class="nf">double_gerr</span><span class="p">(</span><span class="n">xarr</span><span class="p">):</span>
    <span class="k">return</span> <span class="k">lambda</span> <span class="n">p</span><span class="p">:</span><span class="n">xarr</span><span class="o">-</span><span class="n">double_gaussian</span><span class="p">(</span><span class="o">*</span><span class="n">p</span><span class="p">)(</span><span class="o">*</span><span class="n">indices</span><span class="p">(</span><span class="n">xarr</span><span class="o">.</span><span class="n">shape</span><span class="p">))</span></div>
<div class="viewcode-block" id="triple_gerr"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.triple_gerr">[docs]</a><span class="k">def</span> <span class="nf">triple_gerr</span><span class="p">(</span><span class="n">xarr</span><span class="p">):</span>
    <span class="k">return</span> <span class="k">lambda</span> <span class="n">p</span><span class="p">:</span><span class="n">xarr</span><span class="o">-</span><span class="n">triple_gaussian</span><span class="p">(</span><span class="o">*</span><span class="n">p</span><span class="p">)(</span><span class="o">*</span><span class="n">indices</span><span class="p">(</span><span class="n">xarr</span><span class="o">.</span><span class="n">shape</span><span class="p">))</span></div>
<div class="viewcode-block" id="return_param"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.return_param">[docs]</a><span class="k">def</span> <span class="nf">return_param</span><span class="p">(</span><span class="n">xarr</span><span class="p">,</span><span class="n">params</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span><span class="n">negamp</span><span class="o">=</span><span class="bp">False</span><span class="p">):</span>
    <span class="k">if</span> <span class="n">params</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>
        <span class="k">if</span> <span class="n">negamp</span><span class="p">:</span>
            <span class="n">params</span> <span class="o">=</span> <span class="p">[</span><span class="n">xarr</span><span class="o">.</span><span class="n">argmin</span><span class="p">(),</span><span class="mi">5</span><span class="p">,</span><span class="n">xarr</span><span class="o">.</span><span class="n">min</span><span class="p">()]</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">params</span> <span class="o">=</span> <span class="p">[</span><span class="n">xarr</span><span class="o">.</span><span class="n">argmax</span><span class="p">(),</span><span class="mi">5</span><span class="p">,</span><span class="n">xarr</span><span class="o">.</span><span class="n">max</span><span class="p">()]</span>
    <span class="n">pars</span><span class="p">,</span> <span class="n">cov</span><span class="p">,</span> <span class="n">infodict</span><span class="p">,</span> <span class="n">errmsg</span><span class="p">,</span> <span class="n">success</span> <span class="o">=</span> <span class="n">optimize</span><span class="o">.</span><span class="n">leastsq</span><span class="p">(</span><span class="n">gerr</span><span class="p">(</span><span class="n">xarr</span><span class="p">),</span> <span class="n">params</span><span class="p">,</span> <span class="n">full_output</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">pars</span></div>
<div class="viewcode-block" id="return_double_param"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.return_double_param">[docs]</a><span class="k">def</span> <span class="nf">return_double_param</span><span class="p">(</span><span class="n">xarr</span><span class="p">,</span><span class="n">params</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
    <span class="k">if</span> <span class="n">params</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>
        <span class="n">params</span> <span class="o">=</span> <span class="p">[</span><span class="n">xarr</span><span class="o">.</span><span class="n">argmax</span><span class="p">(),</span><span class="n">xarr</span><span class="o">.</span><span class="n">argmax</span><span class="p">()</span><span class="o">+</span><span class="mi">3</span><span class="p">,</span><span class="mf">4.2</span><span class="p">,</span><span class="mf">2.3</span><span class="p">,</span><span class="n">xarr</span><span class="o">.</span><span class="n">max</span><span class="p">(),</span><span class="n">xarr</span><span class="o">.</span><span class="n">max</span><span class="p">()</span><span class="o">/</span><span class="mi">2</span><span class="p">]</span>
    <span class="n">pars</span><span class="p">,</span> <span class="n">cov</span><span class="p">,</span> <span class="n">infodict</span><span class="p">,</span> <span class="n">errmsg</span><span class="p">,</span> <span class="n">success</span> <span class="o">=</span> <span class="n">optimize</span><span class="o">.</span><span class="n">leastsq</span><span class="p">(</span><span class="n">double_gerr</span><span class="p">(</span><span class="n">xarr</span><span class="p">),</span> <span class="n">params</span><span class="p">,</span> <span class="n">full_output</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">pars</span></div>
<div class="viewcode-block" id="return_triple_param"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.return_triple_param">[docs]</a><span class="k">def</span> <span class="nf">return_triple_param</span><span class="p">(</span><span class="n">xarr</span><span class="p">,</span><span class="n">params</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    input parameters: center[1-3],width[1-3],amplitude[1-3]</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">if</span> <span class="n">params</span> <span class="o">==</span> <span class="bp">None</span><span class="p">:</span>
        <span class="n">params</span> <span class="o">=</span> <span class="p">[</span><span class="n">xarr</span><span class="o">.</span><span class="n">argmax</span><span class="p">(),</span><span class="n">xarr</span><span class="o">.</span><span class="n">argmax</span><span class="p">()</span><span class="o">+</span><span class="mi">3</span><span class="p">,</span><span class="n">xarr</span><span class="o">.</span><span class="n">argmax</span><span class="p">(),</span><span class="mf">4.2</span><span class="p">,</span><span class="mf">2.3</span><span class="p">,</span><span class="mi">10</span><span class="p">,</span><span class="n">xarr</span><span class="o">.</span><span class="n">max</span><span class="p">(),</span><span class="n">xarr</span><span class="o">.</span><span class="n">max</span><span class="p">()</span><span class="o">/</span><span class="mf">2.</span><span class="p">,</span><span class="n">xarr</span><span class="o">.</span><span class="n">max</span><span class="p">()</span><span class="o">/</span><span class="mf">5.</span><span class="p">]</span>
    <span class="n">pars</span><span class="p">,</span> <span class="n">cov</span><span class="p">,</span> <span class="n">infodict</span><span class="p">,</span> <span class="n">errmsg</span><span class="p">,</span> <span class="n">success</span> <span class="o">=</span> <span class="n">optimize</span><span class="o">.</span><span class="n">leastsq</span><span class="p">(</span><span class="n">triple_gerr</span><span class="p">(</span><span class="n">xarr</span><span class="p">),</span> <span class="n">params</span><span class="p">,</span> <span class="n">full_output</span><span class="o">=</span><span class="mi">1</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">pars</span>


</div>
<div class="viewcode-block" id="adaptive_collapse_gaussfit"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.adaptive_collapse_gaussfit">[docs]</a><span class="k">def</span> <span class="nf">adaptive_collapse_gaussfit</span><span class="p">(</span><span class="n">cube</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span><span class="n">nsig</span><span class="o">=</span><span class="mi">3</span><span class="p">,</span><span class="n">nrsig</span><span class="o">=</span><span class="mi">4</span><span class="p">,</span><span class="n">prefix</span><span class="o">=</span><span class="s">&#39;interesting&#39;</span><span class="p">,</span>
        <span class="n">vconv</span><span class="o">=</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">x</span><span class="p">,</span><span class="n">xtora</span><span class="o">=</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">x</span><span class="p">,</span><span class="n">ytodec</span><span class="o">=</span><span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">x</span><span class="p">,</span><span class="n">doplot</span><span class="o">=</span><span class="bp">True</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Attempts to fit one or two Gaussians to each spectrum in a data cube and returns the parameters of the fits.</span>
<span class="sd">    Adaptively determines where to fit two Gaussian components based on residuals.  Will fit 3 gaussians if a</span>
<span class="sd">    two-gaussian fit is not better than a certain threshold (specified by nsig), and those fits will be output</span>
<span class="sd">    to images with filename prefix+(coordinate).png.  The 3-gaussian fit parameters will not be returned because</span>
<span class="sd">    the automated fitting is very unlikely to get that part right.</span>

<span class="sd">    inputs:</span>
<span class="sd">    cube - a data cube with two spatial and one spectral dimensions</span>
<span class="sd">    axis - the axis of the spectral dimension</span>
<span class="sd">    nsig - number of sigma over the mean residual to trigger double-gaussian fitting</span>
<span class="sd">           also, cutoff to do any fitting at all</span>
<span class="sd">    prefix - the prefix (including directory name) of the output images from 3-gaussian fitting</span>
<span class="sd">    doplot - option to turn off plotting of triple-gaussian fits</span>

<span class="sd">    vconv,xtora,ytodec - functions to convert the axes from pixel coordinates to ra/dec/velocity coordinates</span>

<span class="sd">    returns:</span>
<span class="sd">    width_arr1,width_arr2,chi2_arr,offset_arr1,offset_arr2,amp_arr1,amp_arr2</span>
<span class="sd">    The Gaussian widths, line centers (in pixel units), amplitudes, and the chi-squared value, not in that order</span>
<span class="sd">    These returns are identical to the returns from double_gaussian, but all components will be zero for the second</span>
<span class="sd">    gaussian in the case of a single-gaussian fit</span>

<span class="sd">    the triple gaussian is guessed to be the double gaussian plus a broad, low-amplitude gaussian.  Ideally this should</span>
<span class="sd">    fit outflows reasonably well, but who knows if it really will.</span>
<span class="sd">    Another option is to fit a negative-amplitude gaussian to account for self-absorption</span>

<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">std_coll</span> <span class="o">=</span> <span class="n">cube</span><span class="o">.</span><span class="n">std</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="n">axis</span><span class="p">)</span>       <span class="c"># standard deviation of each spectrum</span>
<span class="c">#    mad_coll = MAD(cube,axis=axis)</span>
    <span class="n">mean_std</span> <span class="o">=</span> <span class="n">median</span><span class="p">(</span><span class="n">std_coll</span><span class="o">.</span><span class="n">ravel</span><span class="p">())</span>  <span class="c"># median standard deviation (to reject high-signal spectra that have high std)</span>
    <span class="k">if</span> <span class="n">axis</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>                         <span class="c"># force spectral axis to first axis</span>
        <span class="n">cube</span> <span class="o">=</span> <span class="n">cube</span><span class="o">.</span><span class="n">swapaxes</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="n">axis</span><span class="p">)</span>
    <span class="n">width_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>   <span class="c"># define gaussian param arrays</span>
    <span class="n">width_arr1</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>   <span class="c"># define gaussian param arrays</span>
    <span class="n">width_arr2</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>   <span class="c"># define gaussian param arrays</span>
    <span class="n">amp_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>     <span class="c"># define gaussian param arrays</span>
    <span class="n">amp_arr1</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>     <span class="c"># define gaussian param arrays</span>
    <span class="n">amp_arr2</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>     <span class="c"># define gaussian param arrays</span>
    <span class="n">chi2_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>     <span class="c"># define gaussian param arrays</span>
    <span class="n">resid_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>    <span class="c"># define gaussian param arrays</span>
    <span class="n">offset_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>  <span class="c"># define gaussian param arrays </span>
    <span class="n">offset_arr1</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>  <span class="c"># define gaussian param arrays </span>
    <span class="n">offset_arr2</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>  <span class="c"># define gaussian param arrays </span>
    <span class="n">ncarr</span> <span class="o">=</span> <span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="o">*</span><span class="n">nsig</span><span class="p">)</span>        <span class="c"># cutoff: don&#39;t fit no-signal spectra</span>
    <span class="n">starttime</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>              <span class="c"># timing for output</span>
    <span class="k">print</span> <span class="n">cube</span><span class="o">.</span><span class="n">shape</span>
    <span class="k">print</span> <span class="s">&quot;Fitting a total of </span><span class="si">%i</span><span class="s"> spectra with peak signal above </span><span class="si">%f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">ncarr</span><span class="o">.</span><span class="n">sum</span><span class="p">(),</span><span class="n">mean_std</span><span class="o">*</span><span class="n">nsig</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]):</span>      <span class="c"># Loop over all elements for </span>
        <span class="n">t0</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
        <span class="n">nspec</span> <span class="o">=</span> <span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="o">*</span><span class="n">nsig</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>
        <span class="k">print</span> <span class="s">&quot;Working on row </span><span class="si">%d</span><span class="s"> with </span><span class="si">%d</span><span class="s"> spectra to fit&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">nspec</span><span class="p">)</span> <span class="p">,</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">2</span><span class="p">]):</span>
            <span class="k">if</span> <span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="o">*</span><span class="n">nsig</span><span class="p">:</span>
<span class="c">#            if cube[:,i,j].max() &gt; MAD(cube[:,i,j]):            </span>
                <span class="n">pars</span> <span class="o">=</span> <span class="n">return_param</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])</span>
                <span class="n">width_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
                <span class="n">width_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
                <span class="n">amp_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
                <span class="n">amp_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
<span class="c">#                chi2_arr[i,j] = sum(( gerr(cube[:,i,j])(pars) )**2) </span>
                <span class="n">resid_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">gerr</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])(</span><span class="n">pars</span><span class="p">))</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>
                <span class="n">offset_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">offset_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="k">else</span><span class="p">:</span> 
                <span class="n">width_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">chi2_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">resid_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">offset_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
        <span class="n">dt</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span><span class="o">-</span><span class="n">t0</span>
        <span class="k">if</span> <span class="n">nspec</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="k">print</span> <span class="s">&quot;in </span><span class="si">%f</span><span class="s"> seconds (average: </span><span class="si">%f</span><span class="s">)&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">dt</span><span class="p">,</span><span class="n">dt</span><span class="o">/</span><span class="nb">float</span><span class="p">(</span><span class="n">nspec</span><span class="p">))</span>
        <span class="k">else</span><span class="p">:</span> 
            <span class="k">print</span> 
    <span class="n">chi2_arr</span> <span class="o">=</span> <span class="n">resid_arr</span><span class="o">**</span><span class="mi">2</span>
    <span class="n">resids</span> <span class="o">=</span> <span class="n">ma</span><span class="o">.</span><span class="n">masked_where</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">chi2_arr</span><span class="p">),</span><span class="n">chi2_arr</span><span class="p">)</span> <span class="c"># hide bad values</span>
<span class="c">#    residcut = (resids.mean() + (resids.std() * nrsig) )  # Old versino - used standard deviation and mean</span>
    <span class="n">residcut</span> <span class="o">=</span> <span class="p">(</span><span class="n">nanmedian</span><span class="p">(</span><span class="n">chi2_arr</span><span class="o">.</span><span class="n">ravel</span><span class="p">())</span> <span class="o">+</span> <span class="p">(</span><span class="n">MAD</span><span class="p">(</span><span class="n">chi2_arr</span><span class="o">.</span><span class="n">ravel</span><span class="p">())</span> <span class="o">*</span> <span class="n">nrsig</span><span class="p">)</span> <span class="p">)</span> <span class="c"># New version: set cutoff by median + nrsig * MAD</span>
    <span class="n">to_refit</span> <span class="o">=</span> <span class="p">(</span><span class="n">resids</span> <span class="o">&gt;</span> <span class="n">residcut</span><span class="p">)</span><span class="o">.</span><span class="n">astype</span><span class="p">(</span><span class="s">&#39;bool&#39;</span><span class="p">)</span>
<span class="c">#    to_refit[numpy.isnan(to_refit)] = 0</span>
    <span class="n">inds</span> <span class="o">=</span> <span class="n">array</span><span class="p">(</span><span class="n">nonzero</span><span class="p">(</span><span class="n">to_refit</span><span class="p">))</span><span class="o">.</span><span class="n">transpose</span><span class="p">()</span>
    <span class="n">dgc</span><span class="p">,</span><span class="n">tgc</span> <span class="o">=</span> <span class="mi">0</span><span class="p">,</span><span class="mi">0</span>
    <span class="k">print</span> <span class="s">&quot;Refitting a total of </span><span class="si">%i</span><span class="s"> spectra with peak residual above </span><span class="si">%f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">to_refit</span><span class="o">.</span><span class="n">sum</span><span class="p">(),</span><span class="n">residcut</span><span class="p">)</span>
    <span class="n">f</span><span class="o">=</span><span class="nb">open</span><span class="p">(</span><span class="s">&quot;</span><span class="si">%s</span><span class="s">_triples.txt&quot;</span> <span class="o">%</span> <span class="n">prefix</span><span class="p">,</span><span class="s">&#39;w&#39;</span><span class="p">)</span>
<span class="c">#    vconv = lambda x: (x-p3+1)*dv+v0    # convert to velocity frame</span>
    <span class="n">vind</span> <span class="o">=</span> <span class="n">vconv</span><span class="p">(</span><span class="n">arange</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">]))</span>
    <span class="n">xind</span> <span class="o">=</span> <span class="n">arange</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
    <span class="k">for</span> <span class="n">ind</span> <span class="ow">in</span> <span class="n">inds</span><span class="p">:</span>
        <span class="n">i</span><span class="p">,</span><span class="n">j</span> <span class="o">=</span> <span class="n">ind</span>
        <span class="n">doublepars</span> <span class="o">=</span> <span class="n">return_double_param</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])</span>
        <span class="n">old_chi2</span> <span class="o">=</span> <span class="n">chi2_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>
        <span class="n">new_chi2</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">(</span><span class="n">square</span><span class="p">(</span> <span class="n">double_gerr</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])(</span><span class="n">doublepars</span><span class="p">)</span> <span class="p">))</span> 
        <span class="k">if</span> <span class="n">new_chi2</span> <span class="o">&lt;</span> <span class="n">old_chi2</span><span class="p">:</span> <span class="c"># if 2 gaussians is an improvement, use it!</span>
            <span class="n">chi2_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">new_chi2</span>
            <span class="n">width_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
            <span class="n">width_arr2</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span>
            <span class="n">amp_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span>
            <span class="n">amp_arr2</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span>
            <span class="n">offset_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
            <span class="n">offset_arr2</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
            <span class="n">ncarr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">+=</span> <span class="mi">1</span>
        <span class="k">if</span> <span class="n">new_chi2</span> <span class="o">&gt;</span> <span class="n">residcut</span><span class="p">:</span> <span class="c"># Even if double was better, see if a triple might be better yet [but don&#39;t store it in the params arrays!]</span>
            <span class="k">print</span> <span class="o">&gt;&gt;</span><span class="n">f</span><span class="p">,</span><span class="s">&quot;Triple-gaussian fitting at </span><span class="si">%i</span><span class="s">,</span><span class="si">%i</span><span class="s"> (</span><span class="si">%i</span><span class="s">&#39;th double, </span><span class="si">%i</span><span class="s">&#39;th triple)&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">dgc</span><span class="p">,</span><span class="n">tgc</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">tgc</span> <span class="o">%</span> <span class="mi">100</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
                <span class="k">print</span> <span class="s">&quot;Triple-gaussian fitting at </span><span class="si">%i</span><span class="s">,</span><span class="si">%i</span><span class="s"> (</span><span class="si">%i</span><span class="s">&#39;th double, </span><span class="si">%i</span><span class="s">&#39;th triple)&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">dgc</span><span class="p">,</span><span class="n">tgc</span><span class="p">)</span>
            <span class="n">tgc</span> <span class="o">+=</span> <span class="mi">1</span>
            <span class="n">tpguess</span> <span class="o">=</span> <span class="p">[</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">1</span><span class="p">],(</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">+</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span><span class="o">/</span><span class="mf">2.</span><span class="p">,</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">*</span><span class="mf">5.</span><span class="p">,</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span><span class="o">/</span><span class="mf">5.</span><span class="p">]</span>
            <span class="n">triplepars</span> <span class="o">=</span> <span class="n">return_triple_param</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">],</span><span class="n">params</span><span class="o">=</span><span class="n">tpguess</span><span class="p">)</span>
            <span class="n">pars</span> <span class="o">=</span> <span class="p">[</span><span class="n">offset_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">],</span><span class="n">width_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">],</span><span class="n">amp_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]]</span>
            <span class="k">if</span> <span class="n">doplot</span><span class="p">:</span> <span class="c"># if you don&#39;t, there&#39;s really no point in fitting at all...</span>
                <span class="n">ax</span> <span class="o">=</span> <span class="n">axes</span><span class="p">([</span><span class="o">.</span><span class="mo">05</span><span class="p">,</span><span class="o">.</span><span class="mo">05</span><span class="p">,</span><span class="o">.</span><span class="mi">7</span><span class="p">,</span><span class="o">.</span><span class="mi">9</span><span class="p">])</span>
                <span class="n">plot</span><span class="p">(</span><span class="n">vind</span><span class="p">,</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">],</span><span class="n">color</span><span class="o">=</span><span class="s">&#39;black&#39;</span><span class="p">,</span><span class="n">linestyle</span><span class="o">=</span><span class="s">&#39;steps&#39;</span><span class="p">,</span><span class="n">linewidth</span><span class="o">=</span><span class="s">&#39;.5&#39;</span><span class="p">)</span>
                <span class="n">plot</span><span class="p">(</span><span class="n">vind</span><span class="p">,</span><span class="n">gaussian</span><span class="p">(</span><span class="o">*</span><span class="n">pars</span><span class="p">)(</span><span class="n">xind</span><span class="p">),</span><span class="s">&#39;r-.&#39;</span><span class="p">,</span><span class="n">label</span><span class="o">=</span><span class="s">&quot;Single </span><span class="si">%f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span> <span class="p">(</span><span class="n">gerr</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])(</span><span class="n">pars</span><span class="p">))</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span> <span class="p">)</span> <span class="p">)</span>
                <span class="n">plot</span><span class="p">(</span><span class="n">vind</span><span class="p">,</span><span class="n">double_gaussian</span><span class="p">(</span><span class="o">*</span><span class="n">doublepars</span><span class="p">)(</span><span class="n">xind</span><span class="p">),</span><span class="s">&#39;g--&#39;</span><span class="p">,</span><span class="n">label</span><span class="o">=</span><span class="s">&quot;Double </span><span class="si">%f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span> <span class="p">(</span><span class="n">double_gerr</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])(</span><span class="n">doublepars</span><span class="p">))</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span> <span class="p">))</span>
                <span class="n">plot</span><span class="p">(</span><span class="n">vind</span><span class="p">,</span><span class="n">triple_gaussian</span><span class="p">(</span><span class="o">*</span><span class="n">triplepars</span><span class="p">)(</span><span class="n">xind</span><span class="p">),</span><span class="s">&#39;b:&#39;</span><span class="p">,</span><span class="n">label</span><span class="o">=</span><span class="s">&quot;Triple </span><span class="si">%f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span> <span class="p">(</span><span class="n">triple_gerr</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])(</span><span class="n">triplepars</span><span class="p">))</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span> <span class="p">),</span><span class="n">linewidth</span><span class="o">=</span><span class="mi">2</span><span class="p">)</span>
                <span class="n">pars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">vconv</span><span class="p">(</span><span class="n">pars</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
                <span class="n">text</span><span class="p">(</span><span class="mf">1.05</span><span class="p">,</span><span class="o">.</span><span class="mi">8</span><span class="p">,</span><span class="s">&quot;c1 </span><span class="si">%3.2f</span><span class="s"> w1 </span><span class="si">%3.2f</span><span class="s"> a1 </span><span class="si">%3.2f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="nb">tuple</span><span class="p">(</span><span class="n">pars</span><span class="p">),</span><span class="n">transform</span><span class="o">=</span><span class="n">ax</span><span class="o">.</span><span class="n">transAxes</span><span class="p">,</span><span class="n">size</span><span class="o">=</span><span class="s">&#39;smaller&#39;</span><span class="p">)</span>
                <span class="n">dp</span> <span class="o">=</span> <span class="p">[</span> <span class="n">vconv</span><span class="p">(</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="p">,</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span> <span class="n">vconv</span><span class="p">(</span><span class="n">doublepars</span><span class="p">[</span><span class="mi">1</span><span class="p">]),</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span> <span class="n">doublepars</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span> <span class="p">]</span>
                <span class="n">text</span><span class="p">(</span><span class="mf">1.05</span><span class="p">,</span><span class="o">.</span><span class="mi">6</span><span class="p">,</span><span class="s">&quot;c1 </span><span class="si">%3.2f</span><span class="s"> w1 </span><span class="si">%3.2f</span><span class="s"> a1 </span><span class="si">%3.2f</span><span class="se">\n</span><span class="s">c2 </span><span class="si">%3.2f</span><span class="s"> w2 </span><span class="si">%3.2f</span><span class="s"> a2 </span><span class="si">%3.2f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="nb">tuple</span><span class="p">(</span><span class="n">dp</span><span class="p">),</span><span class="n">transform</span><span class="o">=</span><span class="n">ax</span><span class="o">.</span><span class="n">transAxes</span><span class="p">,</span><span class="n">size</span><span class="o">=</span><span class="s">&#39;smaller&#39;</span><span class="p">)</span>
                <span class="n">tp</span> <span class="o">=</span> <span class="p">[</span> <span class="n">vconv</span><span class="p">(</span><span class="n">triplepars</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="p">,</span> <span class="n">triplepars</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span> <span class="n">triplepars</span><span class="p">[</span><span class="mi">6</span><span class="p">],</span> <span class="n">vconv</span><span class="p">(</span><span class="n">triplepars</span><span class="p">[</span><span class="mi">1</span><span class="p">]),</span> <span class="n">triplepars</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span> <span class="n">triplepars</span><span class="p">[</span><span class="mi">7</span><span class="p">],</span> <span class="n">vconv</span><span class="p">(</span><span class="n">triplepars</span><span class="p">[</span><span class="mi">2</span><span class="p">]),</span> <span class="n">triplepars</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span> <span class="n">triplepars</span><span class="p">[</span><span class="mi">8</span><span class="p">]</span>  <span class="p">]</span>
                <span class="n">text</span><span class="p">(</span><span class="mf">1.05</span><span class="p">,</span><span class="o">.</span><span class="mi">4</span><span class="p">,</span><span class="s">&quot;c1 </span><span class="si">%3.2f</span><span class="s"> w1 </span><span class="si">%3.2f</span><span class="s"> a1 </span><span class="si">%3.2f</span><span class="se">\n</span><span class="s">c2 </span><span class="si">%3.2f</span><span class="s"> w2 </span><span class="si">%3.2f</span><span class="s"> a2 </span><span class="si">%3.2f</span><span class="se">\n</span><span class="s">c3 </span><span class="si">%3.2f</span><span class="s"> w3 </span><span class="si">%3.2f</span><span class="s"> a3 </span><span class="si">%3.2f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="nb">tuple</span><span class="p">(</span><span class="n">tp</span><span class="p">),</span><span class="n">transform</span><span class="o">=</span><span class="n">ax</span><span class="o">.</span><span class="n">transAxes</span><span class="p">,</span><span class="n">size</span><span class="o">=</span><span class="s">&#39;smaller&#39;</span><span class="p">)</span>
                <span class="n">title</span><span class="p">(</span><span class="s">&quot;Spectrum at </span><span class="si">%s</span><span class="s"> </span><span class="si">%s</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">ratos</span><span class="p">(</span><span class="n">xtora</span><span class="p">(</span><span class="n">i</span><span class="p">)),</span><span class="n">dectos</span><span class="p">(</span><span class="n">ytodec</span><span class="p">(</span><span class="n">j</span><span class="p">)))</span> <span class="p">)</span> 
                <span class="n">legend</span><span class="p">(</span><span class="n">loc</span><span class="o">=</span><span class="s">&#39;best&#39;</span><span class="p">)</span>
                <span class="n">savefig</span><span class="p">(</span><span class="s">&quot;</span><span class="si">%s</span><span class="s">_</span><span class="si">%s</span><span class="s">.</span><span class="si">%s</span><span class="s">.png&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">prefix</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">))</span>
                <span class="n">clf</span><span class="p">()</span>
            <span class="n">ncarr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">+=</span> <span class="mi">1</span>
            <span class="k">print</span> <span class="o">&gt;&gt;</span><span class="n">f</span><span class="p">,</span><span class="n">triplepars</span>
        <span class="n">dgc</span> <span class="o">+=</span> <span class="mi">1</span>

    <span class="n">f</span><span class="o">.</span><span class="n">close</span><span class="p">()</span>
    <span class="k">print</span> <span class="s">&quot;Total time </span><span class="si">%f</span><span class="s"> seconds for </span><span class="si">%i</span><span class="s"> double and </span><span class="si">%i</span><span class="s"> triple gaussians&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span><span class="o">-</span><span class="n">starttime</span><span class="p">,</span><span class="n">dgc</span><span class="p">,</span><span class="n">tgc</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">width_arr1</span><span class="p">,</span><span class="n">width_arr2</span><span class="p">,</span><span class="n">chi2_arr</span><span class="p">,</span><span class="n">offset_arr1</span><span class="p">,</span><span class="n">offset_arr2</span><span class="p">,</span><span class="n">amp_arr1</span><span class="p">,</span><span class="n">amp_arr2</span><span class="p">,</span><span class="n">ncarr</span>
</div>
<div class="viewcode-block" id="collapse_gaussfit"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.collapse_gaussfit">[docs]</a><span class="k">def</span> <span class="nf">collapse_gaussfit</span><span class="p">(</span><span class="n">cube</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span><span class="n">negamp</span><span class="o">=</span><span class="bp">False</span><span class="p">):</span>
    <span class="n">std_coll</span> <span class="o">=</span> <span class="n">cube</span><span class="o">.</span><span class="n">std</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="n">axis</span><span class="p">)</span>
    <span class="n">mean_std</span> <span class="o">=</span> <span class="n">median</span><span class="p">(</span><span class="n">std_coll</span><span class="o">.</span><span class="n">ravel</span><span class="p">())</span>
    <span class="k">if</span> <span class="n">axis</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">cube</span> <span class="o">=</span> <span class="n">cube</span><span class="o">.</span><span class="n">swapaxes</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="n">axis</span><span class="p">)</span>
    <span class="n">width_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">amp_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">chi2_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">offset_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">starttime</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
    <span class="k">print</span> <span class="n">cube</span><span class="o">.</span><span class="n">shape</span>
    <span class="k">print</span> <span class="s">&quot;Fitting a total of </span><span class="si">%i</span><span class="s"> spectra with peak signal above </span><span class="si">%f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">((</span><span class="n">cube</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">(),</span><span class="n">mean_std</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]):</span>
        <span class="n">t0</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
        <span class="n">nspec</span> <span class="o">=</span> <span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>
        <span class="k">print</span> <span class="s">&quot;Working on row </span><span class="si">%d</span><span class="s"> with </span><span class="si">%d</span><span class="s"> spectra to fit&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">nspec</span><span class="p">)</span> <span class="p">,</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">2</span><span class="p">]):</span>
            <span class="k">if</span> <span class="ow">not</span> <span class="n">negamp</span> <span class="ow">and</span> <span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="p">:</span>
                <span class="n">pars</span> <span class="o">=</span> <span class="n">return_param</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">],</span><span class="n">negamp</span><span class="o">=</span><span class="n">negamp</span><span class="p">)</span>
                <span class="n">width_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
                <span class="n">chi2_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">((</span> <span class="n">gerr</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])(</span><span class="n">pars</span><span class="p">)</span> <span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> 
                <span class="n">offset_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">amp_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
            <span class="k">elif</span> <span class="n">negamp</span> <span class="ow">and</span> <span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span> <span class="o">&lt;</span> <span class="o">-</span><span class="mi">1</span><span class="o">*</span><span class="n">mean_std</span><span class="p">:</span>
                <span class="n">pars</span> <span class="o">=</span> <span class="n">return_param</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">],</span><span class="n">negamp</span><span class="o">=</span><span class="n">negamp</span><span class="p">)</span>
                <span class="n">width_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
                <span class="n">chi2_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">((</span> <span class="n">gerr</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])(</span><span class="n">pars</span><span class="p">)</span> <span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> 
                <span class="n">offset_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">amp_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
            <span class="k">else</span><span class="p">:</span> 
                <span class="n">width_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">chi2_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">offset_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">amp_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
        <span class="n">dt</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span><span class="o">-</span><span class="n">t0</span>
        <span class="k">if</span> <span class="n">nspec</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="k">print</span> <span class="s">&quot;in </span><span class="si">%f</span><span class="s"> seconds (average: </span><span class="si">%f</span><span class="s">)&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">dt</span><span class="p">,</span><span class="n">dt</span><span class="o">/</span><span class="nb">float</span><span class="p">(</span><span class="n">nspec</span><span class="p">))</span>
    <span class="k">print</span> <span class="s">&quot;Total time </span><span class="si">%f</span><span class="s"> seconds&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span><span class="o">-</span><span class="n">starttime</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">width_arr</span><span class="p">,</span><span class="n">offset_arr</span><span class="p">,</span><span class="n">amp_arr</span><span class="p">,</span><span class="n">chi2_arr</span>

<span class="c"># next step: find 2-gaussian fits</span></div>
<div class="viewcode-block" id="collapse_double_gaussfit"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.collapse_double_gaussfit">[docs]</a><span class="k">def</span> <span class="nf">collapse_double_gaussfit</span><span class="p">(</span><span class="n">cube</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">2</span><span class="p">):</span>
    <span class="n">std_coll</span> <span class="o">=</span> <span class="n">cube</span><span class="o">.</span><span class="n">std</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="n">axis</span><span class="p">)</span>
    <span class="n">mean_std</span> <span class="o">=</span> <span class="n">median</span><span class="p">(</span><span class="n">std_coll</span><span class="o">.</span><span class="n">ravel</span><span class="p">())</span>
    <span class="k">if</span> <span class="n">axis</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
        <span class="n">cube</span> <span class="o">=</span> <span class="n">cube</span><span class="o">.</span><span class="n">swapaxes</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="n">axis</span><span class="p">)</span>
    <span class="n">width_arr1</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">width_arr2</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">amp_arr1</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">amp_arr2</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">chi2_arr</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">offset_arr1</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">offset_arr2</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">:])</span>
    <span class="n">starttime</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
    <span class="k">print</span> <span class="n">cube</span><span class="o">.</span><span class="n">shape</span>
    <span class="k">print</span> <span class="s">&quot;Fitting a total of </span><span class="si">%i</span><span class="s"> spectra with peak signal above </span><span class="si">%f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">((</span><span class="n">cube</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">(),</span><span class="n">mean_std</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]):</span>
        <span class="n">t0</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span>
        <span class="n">nspec</span> <span class="o">=</span> <span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span>
        <span class="k">print</span> <span class="s">&quot;Working on row </span><span class="si">%d</span><span class="s"> with </span><span class="si">%d</span><span class="s"> spectra to fit&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">nspec</span><span class="p">)</span> <span class="p">,</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">xrange</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">2</span><span class="p">]):</span>
            <span class="k">if</span> <span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">&gt;</span> <span class="n">mean_std</span><span class="p">:</span>
                <span class="n">pars</span> <span class="o">=</span> <span class="n">return_double_param</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])</span>
                <span class="n">width_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span>
                <span class="n">width_arr2</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span>
                <span class="n">amp_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span>
                <span class="n">amp_arr2</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span>
                <span class="n">chi2_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="nb">sum</span><span class="p">((</span> <span class="n">double_gerr</span><span class="p">(</span><span class="n">cube</span><span class="p">[:,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">])(</span><span class="n">pars</span><span class="p">)</span> <span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> 
                <span class="n">offset_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
                <span class="n">offset_arr2</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">pars</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
            <span class="k">else</span><span class="p">:</span> 
                <span class="n">width_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">width_arr2</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">chi2_arr</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">offset_arr1</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
                <span class="n">offset_arr2</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">nan</span>
        <span class="n">dt</span> <span class="o">=</span> <span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span><span class="o">-</span><span class="n">t0</span>
        <span class="k">if</span> <span class="n">nspec</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="k">print</span> <span class="s">&quot;in </span><span class="si">%f</span><span class="s"> seconds (average: </span><span class="si">%f</span><span class="s">)&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">dt</span><span class="p">,</span><span class="n">dt</span><span class="o">/</span><span class="nb">float</span><span class="p">(</span><span class="n">nspec</span><span class="p">))</span>
    <span class="k">print</span> <span class="s">&quot;Total time </span><span class="si">%f</span><span class="s"> seconds&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">time</span><span class="o">.</span><span class="n">time</span><span class="p">()</span><span class="o">-</span><span class="n">starttime</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">width_arr1</span><span class="p">,</span><span class="n">width_arr2</span><span class="p">,</span><span class="n">chi2_arr</span><span class="p">,</span><span class="n">offset_arr1</span><span class="p">,</span><span class="n">offset_arr2</span><span class="p">,</span><span class="n">amp_arr1</span><span class="p">,</span><span class="n">amp_arr2</span>
</div>
<div class="viewcode-block" id="wrap_collapse_gauss"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.wrap_collapse_gauss">[docs]</a><span class="k">def</span> <span class="nf">wrap_collapse_gauss</span><span class="p">(</span><span class="n">filename</span><span class="p">,</span><span class="n">outprefix</span><span class="p">,</span><span class="n">redo</span><span class="o">=</span><span class="s">&#39;no&#39;</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>

<span class="sd">    redo - if not equal to &#39;no&#39;, then...</span>
<span class="sd">    if collapse_gaussfit succeeded (to the extent that the .pysav files were written),</span>
<span class="sd">    but some part of the file writing or successive procedures failed, re-do those </span>
<span class="sd">    procedures without redoing the whole collapse</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">fitsfile</span> <span class="o">=</span> <span class="n">pyfits</span><span class="o">.</span><span class="n">open</span><span class="p">(</span><span class="n">filename</span><span class="p">)</span>
    <span class="n">dv</span><span class="p">,</span><span class="n">v0</span><span class="p">,</span><span class="n">p3</span> <span class="o">=</span> <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CD3_3&#39;</span><span class="p">],</span><span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CRVAL3&#39;</span><span class="p">],</span><span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CRPIX3&#39;</span><span class="p">]</span>

    <span class="n">cube</span> <span class="o">=</span> <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span>
    <span class="n">cube</span> <span class="o">=</span> <span class="n">where</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">cube</span><span class="p">),</span><span class="mi">0</span><span class="p">,</span><span class="n">cube</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">redo</span><span class="o">==</span><span class="s">&#39;no&#39;</span><span class="p">:</span>
        <span class="n">doubleB</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">collapse_double_gaussfit</span><span class="p">(</span><span class="n">cube</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">))</span>
        <span class="n">doubleB</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">doubleB</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="n">pickle</span><span class="o">.</span><span class="n">dump</span><span class="p">(</span><span class="n">doubleB</span><span class="p">,</span><span class="nb">open</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_doubleB.pysav&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="s">&#39;w&#39;</span><span class="p">))</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">doubleB</span> <span class="o">=</span> <span class="n">pickle</span><span class="o">.</span><span class="n">load</span><span class="p">(</span><span class="nb">open</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_doubleB.pysav&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="s">&#39;r&#39;</span><span class="p">))</span>
    <span class="n">db</span> <span class="o">=</span> <span class="n">doubleB</span>
    <span class="n">gcd</span> <span class="o">=</span> <span class="n">double_gaussian</span><span class="p">(</span><span class="n">db</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">6</span><span class="p">])(</span><span class="n">indices</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">)[</span><span class="mi">0</span><span class="p">])</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">gcd</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_doublegausscube.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="n">gcd</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">gcd</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">doubleResids</span> <span class="o">=</span> <span class="n">cube</span><span class="o">-</span><span class="n">gcd</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">doubleResids</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_doublegaussresids.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>


    <span class="c">#doubleB[4] = (doubleB[4]-v0) / dv + p3-1</span>
    <span class="c">#doubleB[3] = (doubleB[3]-v0) / dv + p3-1</span>
    <span class="n">doubleB</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">doubleB</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">p3</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="o">*</span> <span class="n">dv</span> <span class="o">+</span> <span class="n">v0</span>
    <span class="n">doubleB</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">doubleB</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">p3</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="o">*</span> <span class="n">dv</span> <span class="o">+</span> <span class="n">v0</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">doubleB</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_doublegausspars.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">redo</span><span class="o">==</span><span class="s">&#39;no&#39;</span><span class="p">:</span>
        <span class="n">singleB</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">collapse_gaussfit</span><span class="p">(</span><span class="n">cube</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">))</span>
        <span class="n">pickle</span><span class="o">.</span><span class="n">dump</span><span class="p">(</span><span class="n">singleB</span><span class="p">,</span><span class="nb">open</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_singleB.pysav&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="s">&#39;w&#39;</span><span class="p">))</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">singleB</span> <span class="o">=</span> <span class="n">pickle</span><span class="o">.</span><span class="n">load</span><span class="p">(</span><span class="nb">open</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_singleB.pysav&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="s">&#39;r&#39;</span><span class="p">))</span>
    <span class="n">gc</span> <span class="o">=</span> <span class="n">gaussian</span><span class="p">(</span><span class="n">singleB</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span><span class="n">singleB</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="n">singleB</span><span class="p">[</span><span class="mi">2</span><span class="p">])(</span><span class="n">indices</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">)[</span><span class="mi">0</span><span class="p">])</span>
    <span class="n">singleB</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">singleB</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">p3</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="o">*</span> <span class="n">dv</span> <span class="o">+</span> <span class="n">v0</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">gc</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_singlegausscube.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="n">gc</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">gc</span><span class="p">)]</span><span class="o">=</span><span class="mi">0</span>
    <span class="n">singleResids</span> <span class="o">=</span> <span class="n">cube</span><span class="o">-</span><span class="n">gc</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">singleResids</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_singlegaussresids.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">singleB</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_singlegausspars.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CD3_3&#39;</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CRVAL3&#39;</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CRPIX3&#39;</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CUNIT3&#39;</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CTYPE3&#39;</span><span class="p">)</span>

    <span class="n">doubleResids</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">doubleResids</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">totalDResids</span> <span class="o">=</span> <span class="n">doubleResids</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">totalDResids</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_doublegauss_totalresids.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="n">singleResids</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">singleResids</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">totalSResids</span> <span class="o">=</span> <span class="n">singleResids</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">totalSResids</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_singlegauss_totalresids.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">singleB</span><span class="p">,</span><span class="n">doubleB</span>

</div>
<div class="viewcode-block" id="wrap_collapse_adaptive"><a class="viewcode-back" href="../../agpy.html#agpy.collapse_gaussfit.wrap_collapse_adaptive">[docs]</a><span class="k">def</span> <span class="nf">wrap_collapse_adaptive</span><span class="p">(</span><span class="n">filename</span><span class="p">,</span><span class="n">outprefix</span><span class="p">,</span><span class="n">redo</span><span class="o">=</span><span class="s">&#39;no&#39;</span><span class="p">,</span><span class="n">nsig</span><span class="o">=</span><span class="mi">5</span><span class="p">,</span><span class="n">nrsig</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span><span class="n">doplot</span><span class="o">=</span><span class="bp">True</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    redo - if not equal to &#39;no&#39;, then...</span>
<span class="sd">    if collapse_gaussfit succeeded (to the extent that the .pysav files were written),</span>
<span class="sd">    but some part of the file writing or successive procedures failed, re-do those </span>
<span class="sd">    procedures without redoing the whole collapse</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">fitsfile</span> <span class="o">=</span> <span class="n">pyfits</span><span class="o">.</span><span class="n">open</span><span class="p">(</span><span class="n">filename</span><span class="p">)</span>
    <span class="n">dv</span><span class="p">,</span><span class="n">v0</span><span class="p">,</span><span class="n">p3</span> <span class="o">=</span> <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CD3_3&#39;</span><span class="p">],</span><span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CRVAL3&#39;</span><span class="p">],</span><span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CRPIX3&#39;</span><span class="p">]</span>
    <span class="n">dr</span><span class="p">,</span><span class="n">r0</span><span class="p">,</span><span class="n">p1</span> <span class="o">=</span> <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CD1_1&#39;</span><span class="p">],</span><span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CRVAL1&#39;</span><span class="p">],</span><span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CRPIX1&#39;</span><span class="p">]</span>
    <span class="n">dd</span><span class="p">,</span><span class="n">d0</span><span class="p">,</span><span class="n">p2</span> <span class="o">=</span> <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CD2_2&#39;</span><span class="p">],</span><span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CRVAL2&#39;</span><span class="p">],</span><span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="p">[</span><span class="s">&#39;CRPIX2&#39;</span><span class="p">]</span>
    <span class="n">xtora</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">p1</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">*</span><span class="n">dr</span><span class="o">+</span><span class="n">r0</span>    <span class="c"># convert pixel coordinates to RA/Dec/Velocity</span>
    <span class="n">ytodec</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">y</span><span class="p">:</span> <span class="p">(</span><span class="n">y</span><span class="o">-</span><span class="n">p2</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">*</span><span class="n">dd</span><span class="o">+</span><span class="n">d0</span>
    <span class="n">vconv</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">v</span><span class="p">:</span> <span class="p">(</span><span class="n">v</span><span class="o">-</span><span class="n">p3</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span><span class="o">*</span><span class="n">dv</span><span class="o">+</span><span class="n">v0</span>

    <span class="n">cube</span> <span class="o">=</span> <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span>
    <span class="n">cube</span> <span class="o">=</span> <span class="n">where</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">cube</span><span class="p">),</span><span class="mi">0</span><span class="p">,</span><span class="n">cube</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">redo</span><span class="o">==</span><span class="s">&#39;no&#39;</span><span class="p">:</span>
        <span class="n">adaptB</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">adaptive_collapse_gaussfit</span><span class="p">(</span><span class="n">cube</span><span class="p">,</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span><span class="n">prefix</span><span class="o">=</span><span class="n">outprefix</span><span class="o">+</span><span class="s">&#39;_triple&#39;</span><span class="p">,</span>
            <span class="n">nsig</span><span class="o">=</span><span class="n">nsig</span><span class="p">,</span><span class="n">nrsig</span><span class="o">=</span><span class="n">nrsig</span><span class="p">,</span><span class="n">vconv</span><span class="o">=</span><span class="n">vconv</span><span class="p">,</span><span class="n">xtora</span><span class="o">=</span><span class="n">xtora</span><span class="p">,</span><span class="n">ytodec</span><span class="o">=</span><span class="n">ytodec</span><span class="p">,</span><span class="n">doplot</span><span class="o">=</span><span class="n">doplot</span><span class="p">))</span>
        <span class="n">adaptB</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">adaptB</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="n">pickle</span><span class="o">.</span><span class="n">dump</span><span class="p">(</span><span class="n">adaptB</span><span class="p">,</span><span class="nb">open</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_adaptB.pysav&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="s">&#39;w&#39;</span><span class="p">))</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">adaptB</span> <span class="o">=</span> <span class="n">pickle</span><span class="o">.</span><span class="n">load</span><span class="p">(</span><span class="nb">open</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_adaptB.pysav&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="s">&#39;r&#39;</span><span class="p">))</span>
    <span class="n">db</span> <span class="o">=</span> <span class="n">adaptB</span>
    <span class="n">gcd</span> <span class="o">=</span> <span class="n">double_gaussian</span><span class="p">(</span><span class="n">db</span><span class="p">[</span><span class="mi">3</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">4</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">5</span><span class="p">],</span><span class="n">db</span><span class="p">[</span><span class="mi">6</span><span class="p">])(</span><span class="n">indices</span><span class="p">(</span><span class="n">cube</span><span class="o">.</span><span class="n">shape</span><span class="p">)[</span><span class="mi">0</span><span class="p">])</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">gcd</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_adaptgausscube.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="n">gcd</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">gcd</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">adaptResids</span> <span class="o">=</span> <span class="n">cube</span><span class="o">-</span><span class="n">gcd</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">adaptResids</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_adaptgaussresids.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>


    <span class="c">#adaptB[4] = (adaptB[4]-v0) / dv + p3-1</span>
    <span class="c">#adaptB[3] = (adaptB[3]-v0) / dv + p3-1</span>
    <span class="n">adaptB</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">adaptB</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span><span class="o">-</span><span class="n">p3</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="o">*</span> <span class="n">dv</span> <span class="o">+</span> <span class="n">v0</span>
    <span class="n">adaptB</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">adaptB</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span><span class="o">-</span><span class="n">p3</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span> <span class="o">*</span> <span class="n">dv</span> <span class="o">+</span> <span class="n">v0</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">adaptB</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_adaptgausspars.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CD3_3&#39;</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CRVAL3&#39;</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CRPIX3&#39;</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CUNIT3&#39;</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">header</span><span class="o">.</span><span class="n">__delitem__</span><span class="p">(</span><span class="s">&#39;CTYPE3&#39;</span><span class="p">)</span>

    <span class="n">adaptResids</span><span class="p">[</span><span class="n">numpy</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">adaptResids</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">totalDResids</span> <span class="o">=</span> <span class="n">adaptResids</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>
    <span class="n">fitsfile</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">data</span> <span class="o">=</span> <span class="n">totalDResids</span>
    <span class="n">fitsfile</span><span class="o">.</span><span class="n">writeto</span><span class="p">(</span><span class="s">&#39;</span><span class="si">%s</span><span class="s">_adaptgauss_totalresids.fits&#39;</span> <span class="o">%</span> <span class="n">outprefix</span><span class="p">,</span><span class="n">clobber</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">adaptB</span>
</pre></div></div>

          </div>
        </div>
      </div>
        </div>
        <div class="sidebar">
          <h3>Table Of Contents</h3>
          <ul>
<li class="toctree-l1"><a class="reference internal" href="../../agpy.html">Adam Ginsburg&#8217;s Python Code (agpy)</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../image_tools.html">Image Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../fft_tools.html">AG_fft_tools Package</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../plfit.html">plfit Package</a></li>
</ul>

          <h3 style="margin-top: 1.5em;">Search</h3>
          <form class="search" action="../../search.html" method="get">
            <input type="text" name="q" />
            <input type="submit" value="Go" />
            <input type="hidden" name="check_keywords" value="yes" />
            <input type="hidden" name="area" value="default" />
          </form>
          <p class="searchtip" style="font-size: 90%">
            Enter search terms or a module, class or function name.
          </p>
        </div>
        <div class="clearer"></div>
      </div>
    </div>

    <div class="footer">
      &copy; Copyright 2011, Adam Ginsburg.
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.2pre.
    <script type="text/javascript">

      var _gaq = _gaq || [];
      _gaq.push(['_setAccount', 'UA-6253248-2']);
      _gaq.push(['_trackPageview']);

      (function() {
        var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
        ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
        var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
      })();

    </script>
        
    </div>
  </body>
</html>